Energy Norms and the Stability of the Einstein Evolution Equations 
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The Einstein evolution equations may be written in a variety of equivalent analytical forms, 
but numerical solutions of these different formulations display a wide range of growth rates for 
constraint violations. For symmetric hyperbolic formulations of the equations, an exact expression 
for the growth rate is derived using an energy norm. This expression agrees with the growth rate 
determined by numerical solution of the equations. An approximate method for estimating the 
growth rate is also derived. This estimate can be evaluated algebraically from the initial data, 
and is shown to exhibit qualitatively the same dependence as the numerically-determined rate on 
the parameters that specify the formulation of the equations. This simple rate estimate therefore 
provides a useful tool for finding the most well-behaved forms of the evolution equations. 
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I. INTRODUCTION 

It is well known that the Einstein e< 
written in a variety of forms 

m iMoisji js m rS 

1 g|Tl2§|3ril8riir|30f |32j. In recent 
years a growing body of work has documented the fact 
that these different formulations, while equivalent an- 
alytically, have significantly different stability proper- 
ties when used for unconstrained [|63f numerical evolu- 
tions [|[ 
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, |35j. The most important 
of these differences is the behavior of non-physical solu- 
tions of the evolution equations, which often grow ex- 
ponentially and eventually dominate the desired phys- 
ical solutions. These non-physical solutions could be 
solutions of the evolution equations that violate the 
constraints ( "constraint- violating instabilities") or solu- 
tions that satisfy the constraints but represent some 
ill-behaved coordinate transformation ("gauge instabil- 
ities"). In many cases it is the rapid exponential growth 
of these non-physical solutions, rather than numerical is- 
sues, that appear to be the key factor that limits our 
ability to run numerical simulations of black holes for 
long times [H 



37 1 . For lack of a better term, 



we refer to these non-physical solutions as "instabilities" 
(because they are unstable, i.e., exponentially-growing, 
solutions of the evolution equations), but keep in mind 
that they are neither numerical instabilities nor do they 
represent physics. 

In this paper we explore the use of the energy norm 
(which can be introduced for any symmetric hyperbolic 
form of the evolution equations) to study these instabil- 
ities. We derive an exact expression for the growth rate 
in terms of the energy norm, and verify that the rate de- 
termined in this way agrees with the growth rate of the 
constraint violations determined numerically. We also 
derive an approximate expression for this growth rate 
that can be evaluated algebraically from the initial data 
for the evolution equations. We explore the accuracy 
of this approximation by comparing it with numerically- 
determined growth rates for solutions of a family of sym- 



metric hyperbolic evolution equations. 

In order to compare the analytical expressions for the 
growth rates derived here with the results of numerical 
computations, it is necessary to select some particular 
family of evolution equations with which to make the 
comparisons. Here we focus our attention on the 12- 
parameter family of first-order evolution equations in- 
troduced by Kidder, Scheel, and Teukolsky (KST) [§6j. 
This family of equations has been shown to be strongly 
hyperbolic when certain inequalities are satisfied by the 
12 parameters; however, our expressions for the insta- 
bility growth rates apply only to symmetric hyperbolic 
systems of equations. Therefore we must extend the anal- 
ysis of the KST equations by explicitly constructing the 
symmetrizer (or metric on the space of fields) that makes 
the equations symmetric hyperbolic. We show that such 
a symmetrizer (in fact a four-parameter family of such 
symmetrizers) can be constructed for an open subset of 
the KST equations having only physical characteristic 
speeds. 

We compare numerical evolutions of the symmetric hy- 
perbolic subset of the KST equations with our analyt- 
ical expressions (both exact and approximate) for the 
growth rates of the instabilities. We make these com- 
parisons using two sets of initial data for the evolution 
equations: flat space in Rindler coordinates P8[ , and the 
Schwarzschild geometry in Painleve-Gullstrand coordi- 
nates |S9| [h| pflj . We find that our exact analytic ex- 
pression for the growth rate of the instability agrees with 
the actual growth rate of the constraints in (both fully 
nonlinear and linearized) numerical simulations. This 
agreement provides further evidence that the constraint- 
violating instabilities are real features of the evolution 
equations and not an artifact of using a poor numerical 
algorithm. In addition, the approximate analytical ex- 
pressions for the growth rates derived here are shown to 
have good qualitative agreement with the numerically- 
determined rates. This approximation therefore provides 
a useful tool for finding more well-behaved formulations 
of the equations. Furthermore, the growth rate of the 
instability is shown here to depend in a non-trivial way 
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on the exact "background" solution as well as on the 
particular formulation of the equations. Hence, unfortu- 
nately, it seems likely that it will never be possible to 
find a unique "most stable" form of the equations for the 
evolution of all initial data. 



Here A ka p and F a p may depend on u" but not on Su a . 
We illustrate in Fig. [I] below that the constraint-violating 
instabilities occur in the solutions to these linearized evo- 
lution equations as well as in the solutions to the full 
non-linear equations. 



II. ENERGY NORMS AND RATE ESTIMATES 



A. Energy Evolution 



We limit our study here to formulations of the Einstein 
evolution equations that can be expressed as first-order 
systems 



d t u a + A pduiir = F a 



(2.1) 




Here u a is the collection of dynamical fields, A p and 
F a are (generally complicated) functions of u a , and dt 
and dk are the partial derivatives with respect to the 
time t and the spatial coordinates x respectively. (We 
use Greek indices to label the dynamical fields and Latin 
indices to label spatial co mpo nents of tensors.) Systems 
of equations of the form (2A) are called weakly hyper- 
bolic |42| if nkA ka p has all real eigenvalues for all spa- 
tial one- forms nk , and strongly hyperbolic if in addition 
rikA ka p has a complete set of eigenvectors for all n*,. 
There exists a large literature devoted to a variety of 
representations of the Einstein evolution equations that 
satisfy these conditions | 
H 0, 0, H, H P |a |2J, £4 
In particular, the 12-parameter KST system of equations 
that we use for our numerical comparisons is of this form. 

In order to const ruc t an energy norm, first-order sys- 
tems such as Eq. (2.1) must have an additional struc- 
ture: a "symmetrizer" S a p. First-order systems of evo- 
lution equations are called symmetric hyperbolic [fi2[ (or 
symmetrizable hyperbolic) if there exists a symmetrizer 
which serves as a metric on the space of fields. Such 
a symmetrizer must be symmetric and positive definite 
(i.e. S a pu a u^ > 0, V u a 7^ 0); in addition, it must 
symmetrize the matrices A ka p: S atl A k ^ 1 p = A a o — 
Ah a , V k. In this paper we limit our discussion to 
symmetric hyperbolic formulations. Note that symmet- 
ric hyperbolic systems are automatically strongly hyper- 
bolic, because symmetric matrices n^A k a ^ always have 
real eigenvalues with a complete set of eigenvectors. But 
the converse is not true: strongly hyperbolic systems 
need not be symmetric hyperbolic (except in one spa- 
tial dimension). In Sec. |lll| we construct symmetrizers 
for (an open subset of) the KST equations. 

Let us turn now to the question of the stability of the 
evolution equations. To do this we consider solutions to 
the equations that are close (as defined by the metric 
S a p) to an exact "background" solution u™ Q. Note 
that ti° may be time-dependent. We define 5u a = u a —u" 
to be the deviation of the solution u a from this given 
background solution. The evolution of Su a is determined 
by the linearized evolution equations: 



d t 8u a + A ka pd k 6u f3 = F a p6u . 



For any symmetric hyperbolic system of evolution 
equations, we may define a natural "energy density" and 
"energy-flux" G3, Q associated with 5u a : 



SE = S a pSu a 5u , 
SE k = A a p5u a 5u . 



(2.3) 
(2.4) 



It follows immediately from the linearized evolution equa- 
tions Eq. (2.2) that this energy density evolves as follows, 



d t 6E + V fc SE k = C a p5u a 5u p , 



(2.5) 



where is the spatial covariant derivative associated 
with the (background) three- metric gij, C a p is given by 



C a p — 25 At ( ct F /i ( g) + d t S a p 

+(V9r 1 dk(VgA k a p), 



(2.6) 



and g = det gij is the determinant of the (background) 
spatial metric. Note that C a g , which serves as the source 
(or sink) for the energy in (£Jj), depends on u™ but not 
on Su a . 



B. Exact Expression for the Growth Rate 

Next we explore the possibility of using this energy to 
measure and to estimate the growth rate of instabilities. 
Define the growth rate 1/r of the energy norm to be 



d t \\SE\\ 



r 2\\SE\\ ' 
where the energy norm II^H is defined by 



1 1^1 1 



J 5Eyfg~d 3 x. 



(2.7) 



(2.8) 



Integrating Eq. ( |2.5| ) over a i=constant surface, we obtain 
the following general expression for the growth rate of the 
energy norm: 



1 



2\\SE\\ 



(CapS^Su! 3 -V n 5E n ) y/g~d 3 x. (2.9) 



(2.2) 



We note that Eq. (2.9) is an identity for any solution 
of the equations. The rate 1/r becomes independent of 
time when Su a grows exponentially: 8u a oc e l l r . 

Figure [l] illustrates the equivalence between the energy 
norm measure and the standard measures of the growth 



3 




FIG. 1: Energy norm |j<5-E|j and constraints ||C|| (per unit 
volume) for evolutions of perturbed Schwarzschild initial data 
using three spectral resolutions. Solid curves are \\SE\\ from 
the full non-linear evolution code, and points are from the 
linearized code. Dotted curves are ||C|| from the non- linear 
code. 



rate of the constraint-violating instability. Plotted are 
results from 3D nonlinear numerical evolutions of per- 
turbed Schwarzschild initial data using a particular for- 
mulation of the Einstein evolution equations ]6^ j. The 
solid curves show the evolution of the energy norm \ \5E\ \ 
while the dotted curves show the evolution of the norm 
of the constraint violations. 



||C|| = J {C 2 + C i C i )^gd^x, 



(2.10) 



where C represents the Hamiltonian and Ci the momen- 
tum constraints |36|. The larger points plotted in Fig. |l| 
show the energy norm computed for a numerical solu- 
tion of the linearized evolution equations, indicating that 
the constraint- violating instabilities occur even in the lin- 
earized theory. 

Figure [l] clearly illustrates that the constraint viola- 
tions ||C|| grow at the same rate as the energy norm 
1 1 5 i? 1 1 (until the very end of the simulation when non- 
linear effects become significant). This equality between 
the growth rates is exact for any constraint-violating 
solution having the form 5u a {t,x) = e t ^ T Su a (0,x). 
Our numerical solution approaches this form asymptoti- 
cally. The numerically-determined slopes of these curves, 
1/t(6El) — 0.0489 (from the linear evolution code), 
1/t(8Eml) = 0.0489 (from the non- linear code) and 
1/t(C) = 0.0490 (also from the non-linear code), are also 
in good agreement with the growth rate determined from 
the integral expression in Eq. l/r(f) = 0.0489. 

This agreement shows that the numerical solutions sat- 
isfy the identity in Eq. (2.9). This provides further strong 



than arising from purely numerical problems associated 
with the discrete representation of the solution or the 
time-evolution algorithm. 

The computational domain, boundary conditions, ini- 
tial data, and other details of the numerical evolutions 
sho wn in Fig. |l| are the same as described later in Sec- 



tion IV B. To choose gauge conditions we set the shift 
and the densitized lapse equal to their analytic values 
for all time. Each nonlinear evolution in Fig. |l| is shown 
for three different spectral resolutions, 18x8x15, 24x8x15, 
and 32x8x15 (where the notation N r xNgxN,p represents 
the number of spectral collocation points in the r, 6, and 
4> directions), demonstrating the asymptotic convergence 
of these results. The results for the same three resolu- 
tions using the linear evolution code are indistinguish- 
able from each other in Fig. ^, so only one resolution is 
plotted. These linearized results are also essentially in- 
distinguishable (until very late times) from the highest 
resolution non-linear results. 



C. Approximate Expression for the Growth Rates 

Although Eq. ( |2.9D is an identity, it does not provide a 
particularly useful way to determine 1/r. Its use requires 
the full knowledge of the spatial structure of the unstable 
solution Su a , and this can be determined only by solving 
the equations. Our goal is to obtain a reasonable estimate 
of 1/r without having to solve the evolution equations. 

We first note that if SE k nu > at the boundaries 
(where nu is the outward-directed normal one-form at the 
boundary), one can integrate (2.9) by parts and obtain 



1 1 



C al3 Su a 5u y/gd 3 x. 



(2.11) 



Therefore if A max is the largest eigenvalue A of the equa- 



tion = (C a f3 — XS a p)Su^, then 



- < 
r 



or equivalently, 



\\SE\\<C\\8E {t=0) \\e } 



(2.12) 



(2.13) 



evidence that the constraint-violating instabilities seen 
here are real solutions to the evolution equations, rather 



for some constant C. This argument is often used |44| to 
show that symmetric hyperbolic systems have well-posed 
initial value problems, i.e., that symmetric hyperbolic 
systems have growth rates that are bounded by expo- 
nentials. Because our numerical simulations use bound- 
ary conditions that satisfy SE k nk > (inco ming charac- 
teristic fields are zero), we could use ( 2.12 ) to estimate 
the growth rate. U nfort unately, we find that the upper 
bound provided by ( 2.1 2| ) is typically far greater than the 
actual growth rate, and therefore an estimate based on 
this bound is not very useful. 

Therefore we take a different approach. Without any 
prior knowledge of the structure of the actual unstable so- 
lution, the simplest choice is to assume that the spatial 
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gradients are given approximately by d n 5u a f=a k n du a , 
where k n is a "wavevector" that characterizes the direc- 
tion and magnitude of the gradient of Su a . (Since the ac- 
tual solutions that are responsible for the instabilities in 
the cases we have studied do seem to have a characteristic 
lengthscale, typically the mass of the black hole or some 
other physically distinguished scale, this approximation 
should not be too bad.) In this case the expression in 
Eq. (2.9) for 1/r simplifies to 



1 

T 



1 



211^11 



C a/3 5u a 5u (3 ,/gd 3 x, 



where C a p is given by 



d t S a 







2k n A™ /3 . 



(2.14) 



(2.15) 



Next we limit the 5u a to the subspace of field vec- 
tors that satisfy the boundary conditions. We do this 
formally for these Su a by introducing the projection op- 
erator P a p. This projection annihilates vectors that vi- 
olate the boundary conditions, and leaves those vectors 
that satisfy all the boundary conditio ns un changed p7| . 
Using this projection we re-write Eq. (2.14) as 



1 



2||M5|| 



C a f 3 P a t ,P t3 Ju^6u ,/ ^gd 3 x, (2.16) 



where 5E = S aP P a ^P fj Ju^Su" . We expect that the 
fastest growing solution to the evolution equation will 
be the one driven most str ongly by these "source" terms 
on the right side of Eq. (2.16). Thus we approximate 
(roughly) this most unstable solution as the eigenvector 
Se a of C a pP a ^PP v having the largest eigenvalue: 



(2.17) 



The integrals in Eq. (2.16) are easily evaluated for this 



eigenvector 5e a , giving the following approximate expres- 
sion for the growth rate, 



1 



1 



r 2f5E^gd< 



X, 



(2.18) 



In this approximation then the growth rate of this most 
unstable solution is half the maximum eigenvalue of 
CapP a nP 13 ' v We test the accuracy of this approxima- 
tion in Sec. IV by comparing its predictions with the 



results of numerical solutions to the evolution equations. 



III. EINSTEIN EVOLUTION EQUATIONS 

In this Section we introduce the particular formula- 
tions of the Einstein evolution equations that we study 
numerically to make comparisons with the growth rate 
estimates derived in Sec. 0. The rather general 12- 
parameter family of formulations introduced by Kidder, 
Scheel, and Teukolsky (KST) [^6) is ideal for these pur- 
poses. In this Section we review these formulations and 



derive expressions for the symmetrizer S a p and energy 
density SE (when they exist) that are needed for our 
growth rate estimates. This Section contains a very brief 
review of the derivation and basic properties of the KST 
equations using the notation of this paper, followed by 
a rather more detailed and technical derivation of the 
needed symmetrizer and energy norms. Readers more 
interested in our numerical tests of the growth rate es- 
timates might prefer to skip ahead to that material in 
Sec. Evl. 



A. Summary of the KST Equations 

The KST formulation of the Einstein evolution equa- 
tions begins with the standard ADM Q equations (dis- 
cussed in detail in written as a first-order sys- 
tem for the "fundamental" dynamical variables: Uq — 
{gij, Kij, Dkij}, where gij is the spatial metric, is 
the extrinsic curvature, and Dkij = %dk9ij- We ex- 
press these standard ADM equations in the (somewhat 
abstract) form: 



dt9ij 
d t Kij 



N k d k9lJ + 2g k(i d j) N k 



2NK ih (3.1) 



d t D k 



hi) 



(3.2) 
(3.3) 



where N and N k are the lapse a nd s hift res pec tively. The 
• ■ • on the right sides of Eqs. (3.2) and ( |3.3| ) stand for 
the standard terms that appear in the first-order form of 
the ADM equations, which are given explicitly (up to a 
slight change in notation @) in Eqs. (2.14) and (2.24) of 
KST ^6|. The 12-parameter extension of these equations 
proposed by KST splits naturally into two parts: the first 
part has 5 dynamical parameters (represented by Greek 
letters 7, £, 77, Xj and cr) that influence the dynamics (e.g. 
including the characteristic speeds) of the system in a 
fundamental way, and the second part has 7 kinematical 
parameters (represented by Latin letters z, k, a, b, c, d, 
and e) that merely re-define the dynamical fields. 

The 5-parameter family of dynamical modifications of 
these equations is obtained a) by adding a 4-parameter 
family of multiples of the constraints to the fundamental 
ADM equations, and b) by assuming that a certain den- 
sitized lapse function (which depends on one additional 
parameter) rather than the lapse itself, is a fixed func- 
tion on spacetime. The first modification is obtained by 
addin g m ultiple s of the constraints to the right sides of 
Eqs. (gj) and (|J): 

•■■ + 1 Ng l0 C + CNg ab C a(l3)b , (3.4) 
• ■ ■ + hn N 9k{fij) + iX N 9ijCk, (3.5) 



dtKi 



d t Dk 



kij — 

where 7, £, 77, and x are constants, and the constraints 
(for the vacuum case) are defined by 

C = i[WR- KijRV +K% (3.6) 
d = V,/Y-\ V,A. (3.7) 
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dk9ij " Dkij, 

dkDuj — diDkij. 



(3.8) 
(3.9) 



The second modification conies by assuming that the 
densitized lapse Q, defined by 



Q = log(Ng-°), 



(3.10) 



(rather than the lapse itself) is a fixed function on space- 
time, where g = det(gy) and a is a constant. With these 
modifications the extended ADM equations become a 5- 
parameter family of evolution equations for the funda- 
mental fields Uq. These equations can be written as 



d 



tu 



A k « d k 4 = F«. 



(3.11) 



The quantities A\ a g and Fq are functions of the fields 
Uq and the parameters 7, £, 77, x, and a. We give explicit 
expressions for the Aq"p in Sec. IIIB| and the Appendix. 

The 12-parameter KST family of representations of 
the Einstein evolution equations is completed by adding 
a 7-parameter family of kinematical transformations of 
the dynamical field s to the fundamental representations 
given in Eq. (3.11). These transformations replace 



and Dkij by P« and Mkij according to the expressions: 



p.. — x- 



l kij 

ab w 
z 9i]9 K ab, 



(3.12) 

M kij = [k S a k 6 b S^ + e 8^8% + ag tJ g bc 8 a k + b 9lJ g ab S c k 
+cg k ^ a 3) g hc + dg H ^ } g ab ] D abc . (3.13) 

While these kinematical transformations may seem like 
"trivial" re-paramet erizat ions of the theory, numerical 
results (see e.g. Sec. IV B below) have shown that these 
transformations can have a significant effect on the stabil- 
ity of the equations. The general transformation defined 
in Eqs. (3.12) and (3.13) is a linear transformation from 
the basic dynamical fields Uq to the new set of dynamical 
fields u a = {/,,,. /',,. .1//.,, } : 



(3.14) 



where the transformation T a g depends on the kinemati- 
cal parameters and the metric . The special case with 
k = 1 and z = a = b = c = d = e = corresponds to the 
identity transformation u a = Uq . The explicit represen- 
tations of the inverse transformation (T _1 ) a g are given 
in the Appendix. 

The evolution equations for the new transforme d dy - 
namical fields u a are obtained by multiplying Eq. (3.11) 
by T a 3. The result (after some straightfo rward manipu- 
lation) has the same general form as Eq. (3.11), 



d t u a + A ka 3 d k uP = F a , 



with A ka g and F a given by: 



A 3 



pa = T a gF^ + V aij d t g ij 



2A ka 3 V^D kij , 



(3.15) 



(3.16) 
(3.17) 



where V m i is defined by 

dT a , 



V° 



dgi. 



(T 



(3.18) 



All of the terms on the right side of Eq. ( |3.17j ), except 
the term containing dtgij, depend on the u a (or Uq — 
T~ la gu 13 ) and not its derivatives. In the remaining term, 
the quantity d t gij is to be replaced by the right-hand side 
of Eq. ( |3.l[ ); this introduces the spatial derivatives d k gij, 
which are to be replaced by 2D k ij. Thus in the end F a 
in Eq. ( [3.17 ) is simply a function of u a as required. 

The simple transformation ( 3.16 ) that relates the ma- 
trix Aq^ v of the fundamental representation of the equa- 
tions with A k ^ v ensures that the characteristic speeds of 
the theory are independent of the kinematical parame- 
ters: 

= det (-v8 a B + n k A ka fj ) = det (-v8 a 3 + n k A ka f3 ) . 

(3.19) 

These characteristic speeds (relative to the hypersurface 
orthogonal observers) are also independent of direction 
n k in general relativity. Thus the hyperbolicity of these 
formulations can depend only on the dynamical param- 
eters 7, C, 77, x? an d o. One can show that these char- 
acteristic speeds (relative to the hypersurface orthogonal 
observers) are {0, ±1, ±vi, ±i>2, ±^3}, where 



= 2cr, (3.20) 
= i(7 7 -4? ?C 7-2x-12ax- 37 7 C), (3.21) 
= i(2 + 47-r?-2 7 r; + 2x + 47x-??C)- (3.22) 



In much of the analysis that follows, we will restrict 
attention to the subset of these KST equations where 
the characteristic speeds have only the physical values 



{0, ±1}. This requires v\ 



= «3 = 1 if the theory is 



also to be strongly hyperbolic (see KST |26|]). Thus our 
primary focus will be the 9-parameter family of equations 
in which the parameters <r, 77, and x are fixed by the 
conditions. 



V 

X 



5 + 7C+ IO7 + 67C 

4 + 4C + IO7 + 67C 
~5 + 7C+ IO7 + 67C' 



(3.23) 
(3.24) 

(3.25) 



that are needed to ensure v\ — 1*2 = i>3 = 1- The parame- 
ters 7 and £ are arbitrary so long as 5+7C+IO7+67C 7^ 0. 
Thus the curve 7 = -(7( + 5)/(6C+ 10)is forbidden, but 
7 and £ are otherwise unconstrained |69f . 



B. Symmetric Hyperbolicity of KST 



A first order system such as Eq. ( |3.15 ) is called sym- 
metric hyperbolic if there exists a positive definite sym- 
metrizer S a g such that ^4™^ = S a ^A ni g is symmetric in 
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the indices a and j3, A^p = Ap a , for all field configu- 
rations. We assume here that S a p depends only on the 
spatial metric gij [f70|. ft is convenient to represent the 
symmetrizer as a quadratic form 



dS 2 = S a pdu a du , 



(3.26) 



where du a = {dgij,dPij,dM k ij} denotes the standard 
basis of co- vectors on the space of dynamical fields. The 
most general symmetric quadratic form on the space of 
dynamical fields (which depends only on the metric g^) 
is given by 

dS 2 = A x dG 2 + Bx dP 2 + 2D X dGdP 
+A 2 g lk g jl dg tj dg H + B 2 g lk g° l dP lj dP M 
+2D 2 g lk g^d~g t0 dP kl + Cx 9 kL g la g' b dM {klj) dM {lab) 
+C 2 g kl g la 9 3h [dM k%3 - dM {klj) ] [dM lab - dM (lab) ] 
+C 3 g ij dM}dM 1 j + C A g 13 dM 2 dM 2 
+2C 5 g lj dM]dM 2 . (3.27) 

Here dG and dP are the traces of dgij and dPij respec- 
tively, and dgij and dPij are their trace- free parts. The 
two traces of dM k ij are defined by 

dM\ = dM ijk g jk (3.28) 
dMj = dM kij gi k , (3.29) 

and its trace-free part, dM k ij, is 

dM klJ = dM kij + ±[dM\ igj)k -2dM\g l3 

+dMlg ij -SdM\ i g j)h ]. (3.30) 



This quadratic form, Eq. (3.27), is positive definite iff 
{Ax,A 2 ,Bx,B 2 ,Cx,C 2 ,C3,C4} are all positive and Cf < 
C 3 C 4 , D 2 < AxBx, and D 2 < A 2 B 2 . (The signs of C 5 , 
Dx and D 2 are irrelevant.) 

The question now is whether the constants A a, Ba, 
G a and Da can be chosen to make the A^a symmetric 
in a and [3. In order to answer this we need explicit 
expressions for the matrices A na p. Quite generally these 
matrices are determined for these equations by a set of 
12 constants fiA and va, which in turn are determined 
by the 12 parameters of the KST formulations. (We give 
the explicit expressions for fiA and va in terms of the 
KST parameters in the Appendix.) The equations that 
define the A na p in terms of these constants are: 



9t9ij 
dtP tJ 



N n d n9ij , 
N n d n K l3 



(3.31) 



N 



kt j 



+^g bc 8 n {l 8 d j) + h i i g cd 8 n {l 8 b l) 
•/';,//""'//'".'/,, + ^g nb g cd g l3 ] d k M bcd , (3.32) 
N n d n M kij -N[ Vl 8 n k 8\8 c 3 + v 2 6 n {l S» j} 8 c k 
+^g nb g k ( l S c 3) + v 4 g nb gij S c k 



+^g bc 9 k ( i 8 n j) + v 6 g bc gij 5 n k d n P bc , (3.33) 



where ~ means that only the principal parts of the 
equations have been represented explicitly (dtu a + 
A na pd k u> 3 ~0). 

We now evaluate A n a& = S afl A n »p and A n pa = 
Sn u A nfJ ' r y using the expressions in Eqs. ( Ujj ) through 
( |3.33|) . After lengthy algebraic manipulations, we find 
that A™ a is symmetric iff Dx = D 2 = and the following 
constraints are satisfied by the constants Ba and C a- 



= 



B 2 (fix + ^)-Cx(vx+v 2 ), (3.34) 

B 2 {2fix-fi 2 )-C 2 {2vx-V2), (3.35) 

Bx(3/j,x + 3^4 + 9/i 6 ) 

— Cs(3ux + v 2 + U3 + 3m + 3^5 + 9m) 

-C 5 (m + 2v 2 + 2m + va + 6v 5 + 3m), (3.36) 

Bx{2>^ 2 + 3^3 + 9/x 5 ) 

— C 5 (3m + m + m + 3m + 3^5 + 9m) 

-C 4 (m + 2v 2 + 2m + v 4 + 6m + 3m), (3.37) 

B 2 (-2fj,x + 3^2 + 10/i 4 ) - C 3 (10u 2 + 10m + 30m) 

-C 5 (10m + 5v 2 +20m + 10m), (3.38) 

b 2 ( 6^1+ /i 2 + i0M3)-C 5 (i0m + i0m + 30m) 

-C 4 (10m + 5m + 20m + 10m)- (3.39) 



This is a system of six linear equations for the seven 
parameters Ba and C a- Thus we expect there to exist 
solution(s) to these equations for almost all choices of 
the ha and va- However, there is no guarantee that such 
solutions will satisfy the positivity requirements needed 
to ensure that S a p is positive definite. 

We divide into two parts the question of determining 
when solutions to Eqs. (3.34) through (3.39) exist that 
satisfy the appropriate positivity conditions: first the 
question of when a positive definite symmetrizer, S^p, 
exists for the subset of the KST equations whose dy- 
namical fields are the fundamental fields Uq , and second 
the question of when this fundamental symmetrizer can 
be extended to a symmetrizer for the full 12-parameter 
set of KST equations. We consider the second ques- 
tion first. Assume that for a given set of dynamical pa- 
rameters there exists a positive definite S® g such that 



S^A^p = S pi A^ a . Now define S af: 



Q rp—lfj, q0 rp — 1 V 



(3.40) 



One can verify, using Eq. ( 3.16 ), that this S a p sym- 



metrizes A na p. Further it follows, using Eq. (3.14), that 
this S a p is positive definite, 



S^U 13 = SSfl«0«0 > °- 



(3.41) 



since S^p is assumed to be positive definite. In the Ap- 
pendix we give explicit expressions for the constants Ba 
and Ca that define S a p in terms of the constants B A and 
C A that define S^g and the parameters that define the 
transformation T~ la p. 

Now we return to the first, and more difficult, question: 
when does there exist a positive definite symmetrizer S^p 
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for the subset of the KST equations whose dynamical 
fields are the fundamental fields Uq ? At the present time 
we have not solved this problem completely. Rather, we 
restrict our attention to an interesting (perhaps the most 
interesting) subset of these KST equations in which the 
characteristic speeds vi, v 2 and V3 are all the speed of 
light: v\ = t>2 = U3 = 1. The restrictions that these 
conditions place on the dynamical parameters are given 
in Eqs. ( [3.23 ) through ( [3.25| ). Each of these systems is 
strongly hyperbolic. 

One can now evaluate the ha and va a ppro priate for 
t his s ubset of KST equ ation s using E qs. (A. 5 ) through 
( |A.12| ) along with E qs. (|3.23| ) throu gh (|3.25| ). Substitut- 
ing these into Eqs. (3.34) through (3.39) gives the sym- 
metrization conditions for these equations. These con- 
ditions are degenerate in this case, reducing to only five 
independent equations. Solving these five symmetriza- 
tion equations for the C\ in terms of the B A gives 



C? = 



C° 3 



c° 



H3 + C)£ 2 °, 



9(l+7) 2 5? 



= (2 + 3 7 ) 2 B 1 



3(C~5) 2 
10(5 + 3C) 
(9C - 5) 2 



B 2i 



Bo, 



30(5 + 3C) 1 
-3(1 +7) (2 + 37)5? 

(C-5)(9C-5) Q 
10(5 + 30 2 ' 



(3.42) 
(3.43) 

(3.44) 
(3.45) 



These equations guarantee that {C?, C 2 , C3 , C4 } are pos- 
itive for any positive and B 2 if and only if 



> C > 



(3.47) 



The only remaining condition needed to establish sym- 
metric hype rbolic ity is to e nsure that C3C4 — (C5) 2 > 0. 



Using Eqs. fl3.44|) through (|3.46|) it follows that 



(C° 5 f = 



3(5 + 7C + IO7 + 67C) ; 
10(5 + 30 



B^Bl. (3.48) 



Thus the right side is positive whenever B® and B 2 are 
positive iff C > ~f and 5 + 7( + IO7 + 67C + 0. The first 
of these c ondi tions along with Eq. ( 3.42 ) demonstrates 
that Eq. ( |3.47[ ) is the necessary and sufficient constraint 
on the parameters {0 7} to ensure symmetric hyperbol- 
icity. The second of these conditions was als o req uired 
to ensure that the parameters r\ and x m Eqs. (3.24) and 
(3.25) are finite, so it does not represent a new restriction. 

Thus a large open set of this two-parameter family 
of the fundamental KST representations of the Einstein 
evolution equations is symmetric hyperbolic. And per- 
haps even more surprising, the complimentary subset of 
these strongly hyperbolic equations (i.e. when £ > 
or £ < — |) is not symmetric hyperbolic ]7l| . Fu rther , 
the extension of this two-parameter family via Eq. ( 3.40| ) 
produces a nine-parameter family of strongly hyperbolic 



representations the Einstein equations. A large open sub- 
set of this nine-parameter family is symmetric hyperbolic 
(i.e. those that are extensions of the symmetric hyper- 
bolic fundamental representations), while its compliment 
is not symmetric hyperbolic. 

The construction used here to build a symmetrizer 
S a g for the KST equations has succeeded unexpectedly 
well. We found not just a single symmetrizer, but in 
fact a four-parameter family of such sym metri zers. Us- 
ing the expressions for the Ca from Eqs. (3.42) through 
( |3.46D , we see that the symmetrizer S°p is a sum of 
terms that depend linearly on each of the four parameters 
{Ai, A 2l Bi , B 2 }. Thus we may write this symmetrizer 
in the form: 



Q° — A QlO 1 A c20 1 R 0c30 



B° 2 Si% (3.49) 



where the "sub-symmetrizers" S^g represent a set of pos- 
itive definite (under the conditions needed for symmet- 
ric hyperbolicity established above) metrics on mutually 
orthogonal subspaces in the space of fields u a . The ex- 
pressions for these sub-symmetrizers are 



S^ g du a du^ 
S™ g du a du p 



g l3 g kl dgi 3 dg k h 



(3.50) 

ijuyki, (3.51) 

Sl°pdu a du f:) = dP 2 + g lj [3(1 + j)dM} - (2 + 3 7 )rfM 2 ] 
x [3(1 + i)dM] - (2 + fry)dMf] , (3.52) 



(3.46) S% du a du p = g ik gi l dP l3 dP kl ~ Qg k ' 



g la gJ»dM ikij) dM {lab) 



+ [30(5 + 30fV [3(C ~ 5)dM/ - (9C - 5)dM 2 ] 
x [3(C - 5)dMj - (9C - 5)dMf] 

+i(3 + 09 hl g ia g jb [dM kij - dM (Hj) ] 



: [dMiab - dM { 



(lab)] 



(3.53) 



The dimensions of the corresponding sub-spaces are 
{1,5,4,20}. Just as the fundamental symmetrizer S^ g 
is related to the more general sym metrizer S a g by the 
transformation given in Eq. ( 3.40| ), so the fundamen- 
tal sub-symmetrizers are related to the general sub- 
symmetrizers by, 



D aj3 



rp — 1 nAOrp— 1 V 
1 a^fiu 1 P- 



(3.54) 



We note that this transformation leaves the first two sub- 
symmetrizers unchanged: S^ g — S^° g and S 2 g = S^ g . 

One interesting subset of the nine-parameter family of 
KST equations studied here is the two-parameter "gener- 
alized Einstein-Christoffel" system studied extensively by 
KST (^6). In the language used here this two-parameter 
family is defined by 



c 

d 
b 
c 
e 



-k = -1, 
7 + 3z + 3z7, 
—7 — 2z — 37Z, 
-d = 2, 
0. 



(3.55) 
(3.56) 
(3.57) 
(3.58) 
(3.59) 



8 



Substitut ing t hese param eter value s into Eqs. (3.42) 
through fl3.46| ) and ( |A.38[ ) through ( |A.44[ ), we find the 
constants that determine the symmetrizer to be 



C 3 



{l + iz) 2 B 1 , 

C\ = C-2 = B 2l 

Bi + -kB 2l 
~iB 2 , 



(3.60) 
(3.61) 
(3.62) 
(3.63) 
(3.64) 



where B\ and B 2 are arbitrary positive constants. For 
this case the sub-symmetrizers S^g and S*g have the 
simple forms: 



Sl p du a du f3 = dP 2 + g ij dMldMj, (3.65) 

S A a!3 du a du = g kl g la g lb dM kl] dM lab - \g^dM]dM) 

1 

ij 1 



+g lk g 3l dP i3 dP k i - idP 2 . (3.66) 



We also point out that the symmetrizer becomes partic- 
ularly simple in this generalized Einstein-Christoffel case 
by taking A\ = B\ = | and A 2 = B 2 = 1: 

dS 2 = S a gdu a du p = g lk g ll d gi] dg kl +g lk gi l dP l3 dP kl 
+g kl g la g 3b dM kl] dM lab . (3.67) 



present we have not found a use for this unexpected abun- 
dance of symmetrizers and energy norms. In our numeri- 
cal work below we choose (fairly arbitrarily) one member 
of this family to compute our growth-rate estimates. 

In our earlier discussion we found it useful to analyze 
the symmetrizers associated with these equations in two 
steps: first, to consider the symmetrizers associated with 
the fundamental representations of the theory and sec- 
ond to work out how the general symmetrizer can be 
obtained from the fundamental representation by per- 
forming a suitable transformation. Here we find it useful 
to consider the corresponding questions for the energies. 
First we recall that the dynamical fields u a are related 
to the fundamental fields by the transformation, 



rpCt 



(3.72) 



where the matrix T a g [defined in Eqs. fl3.12| ) and ( |3.13[ )] 
depends on the kinematical parameters and the metric 
gij . Thus perturbations Su a are related to perturbations 
in the fundamental fields by: 



(3.73) 



where 



This represents a kind of "Euclidean" metric on the space 
of fields, in which the symmetrizer is just the sum of 
squares of the components of the dynamical fields. 



C. The KST Energy Norms 

The symmetrizers for the KST equations may be writ- 
ten as arbitrary positive l inear combinations of the sub- 
symmetrizers as in Eq. ( 3.49 ). Since the equation for 
C a g, and hence the equation for the evolution of the en- 
ergy is linear in S a g, it follows that there are in fact four 
independent energy "sub-norms" for the KST equations. 
Each is defined using the corresponding sub-symmetrizer: 



5E A 



S^ p 5u a 5u p , 



8E\ = A n » B 5u a 5u 



(3.68) 
(3.69) 



It follows that these energies each satisfy evolution equa- 
tions analogous to Eq. (Ej|): 



d t SE A + V n 5E n A = C^5u a Su^ 



(3.70) 



where C£p is given by 



a 



a/3 



dtS^B 



+(v?r 1 3n(vs^^v). (3.7i) 



Each of these energies gives the same growth rate r 
from Eq. (2.9) for solutions that grow exponentially. At 



v a g = r- lQ 7 a /3 T 7 e 



(3.74) 



One can now work out the tra nsform ation prop erties 
of the energy and flux using Eqs. ( 3 . 1 6| ) and ( 3.40j ): 



6E = S a g5u a 8u p , 

= SE + S° a ,V$5u%8ug, 



where is defined as 



(3.75) 



(3.76) 



Thus the energy 5E is just the fundamental energy SEq — 
S^SuqSuq plus terms proportional to V£g. A similar 
argument leads to the transformation for the energy flux: 

SE n = SE£ + A n Qtlv V^5u^84. (3.77) 

We note that the only dependence of the energy and flux 
on the kinematical parameters comes through V£a and 
hence V^ a . 

Finally, we note that the expression for the transfor- 
mation of the terms C a gSu a 5u^ can be obtained by a 
similar calculation. The result is 



C a g8u a 5u 



°a/3 



HVar^n^A^V^+Eag 



SuqSuq, 
(3.78) 



9 



where the term E a f3 is defined by 

E a/3 5u%6u$ = -2 (5» + V\) 5u«S%V\A™ f3 d n 54 
+2 (*£ + V» a ) Su^A^V^dJu^ (3.79) 

This expression is ob tained by stra ightforward calcula- 
tion using E qs. (3.75 ) through ( 3.77 ). Note that the left 
side of Eq. (3.79) is a quadratic form in Suq while the 
right side depends on the derivatives d n 6v,Q. To under- 
stand this we use the fact that V^ a is non-zero only when 
the index a corresponds to one of the spatial metric com- 
ponents^. e. V a g5u n = V m i5gij. This follows from 



Eqs. ( p.l8|) and ( |3.74| ) and the fact that T a p depends on 
<7ij but none of the other dynamical fields. It follows that 
the term V v a A™ pdnSug includes only the derivatives of 
fields that are present in the me tric e volution equation. 
But this evolution equation, Eq. (3.31), includes only the 
derivatives of the metric along the shift vector, and so it 
follows that V u a A^ a p = -N k V v p. Thus the expression 
for E a may be written in the form: 



E a/3 5u%6u% 



4(6^ + V^ a )(A k ^ + N k S^) 5u% 
xV vij 5D kij (3.80) 



The last terms in this expression came from the term 



VpdkSuQ, which [using Eq. ( 3.18; )] depends only on 



the spatial derivatives of the metric perturbations Sgij : 
V v pd k 8v% = \ ' '' <)),')</,, = 2V vij 6D kij . Thus the right 
side of Eq. ( 3.8C ) is a quadratic form in the Suq as re- 
quired. Finally we note that C a p8u a 8u^ like the energy 
and energy-flux, depends on the kinematical parameters 
only through V a p. The expression for the transforma- 
tion of the matrix C a p that appears in our approximate 
expressions for the insta bility growth-rates, Eq. (2.16), 
follows directly from Eq. (3.7£): 



C a f35u a 5w 



+E a0 -2k n A^ v V% 



6u%5u%. (3.81) 



IV. NUMERICAL TESTS 

In this section we compare the approximate expres- 
sions for the instability growth rates developed in Sec. Q 
with growth rates determined directly from numerical 
solutions of the Einstein evolution equations. For this 
study we use the 2-parameter subset of the KST equa- 
tions, discussed in Sec. 



hib 



referred to as the gen- 
eralized Einstein-Christoffel system p6fl . Since we do 
not yet understand the mean ing of the different energy 
norms developed in Sec. [II C, we limit our consideration 



here to the norm computed from the symmetrizer with 
Ax = = 1/3 and A 2 = B$ = 1. This choice is the 
closest analog we have of the simple "Euclidean" metric 
of Eq. ( 3. 67] ) for these systems. 



We examine the accuracy of the approximate expres- 
sion for the growth rate, Eq. ( 2. IS ), by examining the 
evolution of perturbations about two rather different 
background spacetimes: flat spacetime in Rindler coordi- 
nates J3^], and the Schwarzschild geometry in Painleve- 
Gullstrand coordinates J3^, fl0| ^T|. In each of these 
cases the full 3D numerical evolution of these equations 
has constraint-violating (and possibly gauge) instabili- 
ties, and the approximate expressions for the growth 
rates using our new formalism are simple enough that 
we can evaluate them in a straightforward manner (even 
analytically in some cases) . Thus we are able to compare 
the estimates of the growth rates with the full numerical 
evolutions in a systematic way. 



A. Rindler Spacetime 

Flat spacetime can be expressed in Rindler coordinates 
as follows, 



ds 2 = -x 2 dt 2 



dx + dy + dz 2 



(4.1) 



In this geometry the dynamical fields have the following 
simple forms, 



9ij = 

Kii = 

Dkij = 

N = 

N k = 



(4.2) 
(4.3) 
(4.4) 
(4.5) 
(4.6) 



where Sij represents the Euclidean metric in Cartesian 
coordinates. Figure || illustrates that even this simple 
representation of flat space is subject to the constraint 
violating instabilities. This figure shows the evolution of 
the norm s of each of the constraints defined in Eqs. ( |3.6| ) 
through ( |3.9[ ). This figure also illustrates that all of the 
constraints grow at the same exponential rate in our sim- 
ulations. 



The simplicity of the expressions (4.2) through (4.6) al- 
lows us to evaluate the various quantities needed to make 
our stability estimates in a reasonably straightforward 
way. The first consequence of this simple form is that 
the unperturbed fundamental fields consist only of met- 
ric fields: Uq = {5^,0,0}. This fact considerably simpli- 
fies many of the needed expressions. In particular, the 
quantity daT^ v Uq vanishes identically for this geometry 
because T M „, when v has values that correspond to the 
components of the metric, is just the identity transforma- 
tion and hence independent of the dynamical fields. Thus 
the matrix V^ a defined in Eq. ( |3.74 ) vanishes identi- 
cally for the Rindler geometry. It follows that for Rindler 
SE = SE , 5E n = 6E$, and C^S^Sul 3 = C^5u%Sv%. 
Thus the timescale 1/r associated with the instability 
in the Rindler geometry is completely independent of 
the kinematical parameters of the representation of the 
evolution equations. This independence of 1/r on the 
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FIG. 2: Solid curve shows the evolution of the sum of the inte- 
gral norms of all of the constraints. Dotted curves show the in- 
dividual contributions from the various constraints: CkiijC k l \ 
C 2 , CkijC h%i , and CkC k (in that order from largest to smallest). 



FIG. 3: Solid curve shows (half) the largest eigenvalue of C a f) 
for the Rindler spacetime (with k n — 0). Dots give values for 
the actual growth rate of the short timescale instability as 
determined by numerical solution of the linearized evolution 
equations with f = — 1. 



kinematical parameters has been verified numerically for 
the 2-parameter generalized Einstein-Christoffel subset 
of the linearized KST equations. 

Next we evaluate the simple estimate of the growth 
rate of the fastest-growing mode using Eq. (ElJ). First 
we evaluate the quadratic form C^^Suq Su^. For the 
Rindler geometry with k n — this quantity has a rea- 
sonably simple form: 

C%5u a 8u fi = 

28K ij [(C + 3)x k 8D ik j - 2x k SD klj - 2xS gij ] 

+x k 8K\{{Q 1 - 2)«v - [Qi + 2 (C - 7)]«V} 

+x k 6K\ [(Qi + 2)SD jk j - Q 2 «V] , (4-7) 

where and Qi and Q 2 are given by 



Qx = 2 7 (5 + 3 7 ) + 10 
Q 2 = 6 7 (2 + 7 ) + 8 



1-C 

3C + 5' 
C + 5 



3C + 5 



(4.8) 
(4.9) 



This symmetric quadratic form is simple enough that its 
eigenvalues and eigenvectors can be determined analyt- 
ically j72|. All of these eigenvalues depend on the two 
dynamical parameters {C,7} but not on any of the kine- 
matical parameters. The maximum eigenvalue also de- 
pends on position in Rindler space, according to A max = 
/(Cj7) + 4x 2 , where x is the Rindler coordinate. This 
eigenvalue has the (approximate) minimum value = 
10.198 + 4x 2 at the point {C7} = {-0.135,-1.382}. In 
our growth rate estimates we evaluate the eigenvalue at 
the point in our computational domain where A max has 
its maximum value, which in our case is at x = 1. 

We illustrate in Fig. [|the dependence of 5A max on the 
parameter 7 (for C, — —1). Our simple analysis predicts 



that the best-behaved form of the evolution equations 
should be obtained for parameter values near this min- 
imum. We also plot in Fig. || numerically-determined 
points that represent the growth rates of the instabil- 
ity. These points were evaluated from the growth rate 
of the energy norm for evolutions of the linearized KST 
equations using a ID pseudospectral method on the do- 
main x S [0.1,1]. The numerical method (but not the 
system of equations) is identical to the one described 
in except here our spatial coordinate is the Carte- 
sian coordinate x rather than the spherical coordinate 
r. All evolutions are performed at multiple resolutions 
(see Fig. |J) to test convergence, and the growth rates 
that we quote in the text and in the figures are taken 
from the converged solution. At the boundaries we de- 
mand that incoming characteristic fields have zero time 
derivative, and we do not constrain the outgoing and 
nonpropagating characteristic fields. The initial data is 
taken to be the exact solution plus a perturbation of the 
form yl e -(*-0-55) /to t k at j g added ^ eacn f ^ ne dy- 
namical variables (gij,Kij,M k ij). The perturbation am- 
plitude A for each variable is a randomly chosen number 
in the interval (— 10~ 8 , 10~ 8 ). Note that this perturba- 
tion violates the constraints. The gauge fields (3 % and Q 
are not perturbed. We measure the growth rates of the 
energy norm during the very earliest parts of the evolu- 
tions, before the initial Gaussian pulse can propagate to 
the edge of the computational domain. To do this we 
use an extremely narrow pulse w = 0.0125, so that there 
is sufficient time to measure the growth rate before even 
the tail of the pulse reaches the boundary (one could also 
widen the computational domain, but this is equivalent 
to changing the pulse width because the Rindler solution 
is scale-invariant). 
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Figure |^ illustrates that the analytical estimate con- 
sisting of 1/r « |A max gives a reasonably good estimate 
of this initial growth rate of the instability. And the loca- 
tion of the optimal value of the parameter 7, where the 
growth rate of the instability is minimum, also agrees 
fairly well with the location of the minimum of iA max . 

However, the short timescale instabilities illustrated 
in Fig. H are not our primary concern here. Figure |] 
illustrates the evolution of the energy norm for one of 
the evolutions discussed above (with 7 = |, C = — F 
w = 0.025). This shows that once the unstable initial 
pulse crosses and leaves the computational domain (in 
about half a light crossing time, or t ~ 0.6), the solution 
grows much less rapidly. While the short term instabil- 
ity does not seriously effect the long term stability of the 
code (unless it is large enough that the code blows up in 
less than a crossing time), Fig. ^ illustrates that there are 
often other instabilities that grow more slowly, but which 
can contaminate and eventually destroy any attempts to 
integrate the equations for long periods of time. Figure |] 
also illustrates the equivalence between evolving the sys- 
tem using a 1+1 dimensional finite difference code and 
a pseudo-spectral code. The finite difference code uses a 
first-order upwind method (see, e.g. ]4(|) in which the 
fundamental variables are decomposed into characteristic 
fields. Both the finite difference and the pseudo-spectral 
methods use the same initial data, boundary conditions, 
and gauge conditions. Three different spatial resolutions 
are illustrated for each method, and the highest resolu- 
tion runs essentially coincide. This agreement provides 
additional evidence that the instabilities discussed here 
are features of the analytical evolution equations, and are 
not numerical in origin. 

In our estimates, we will now attempt to filter out 
the less interesting short term instabilities by imposing 
boundar y con ditions on the trial eigenfunctions 5u a used 
in Eq. (2.14). These boundary conditions are imple- 
mented using the projection operator P^ v in Eq. (2.16). 
For the case of the Rindler geometry this projection op- 
erator is constructed to annihilate both the ingoing and 
outgoing propagati ng mo des. Thus the growth rate esti- 
mate given in Eq. ( 2.18| ) is implemented by finding the 
largest eigenvalue of C a pP a ^P 13 „ that is projected onto 
the subspace of non-propagating modes (i.e. the eigen- 
vectors of ikA ka p having eigenvalue zero, as measured 
by a hypersurface orthogonal observer). We find that 
the largest eigenvalue of this C a pP a ^P^ v is zero for all 
values of {C>7}: for all values of the wavevector k n , and 
for all values of x: the non-propagating modes of Rindler 
are all stable according to this estimate. 

This estimate of the long-term growth rate in Rindler 
is shown as the solid curve in Fig. ||[ The points in 
Fig. ^| represent growth rates determined numerically 
over many many light crossing times (for these evolu- 
tions, we have no need for an extremely narrow pulse; we 
use w = 0.1 so that we can run at lower resolutions). We 
observed that these equations are in fact stable for most 
values of the parameter 7 over these timescales. The only 
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FIG. 4: Evolution of the energy norm for perturbations of 
Rindler space. The solid curves show the evolutions based on 
a three different resolutions (101, 401 and 6401 grid points) 
of a finite difference version of the code, while dotted curves 
show the same evolution using three different spectral resolu- 
tions (64, 96, and 128 basis functions). A magnified view of 
the first part of the evolution (upper left) illustrates the short 
timescale instability whose growth rate is shown in Fig. H. 




FIG. 5: Solid curve illustrates the simple analytical growth 
rate estimate 1/r = for the Rindler spacetime. Points 
are growth rates determined numerically for long-term (many 
light crossing times) evolutions with f = — 1. Dotted curve 
represents the short-term instability growth rate estimate. 



long term instabilities that we observed in Rindler occur 
for values of 7 near 0.5, where the short-term instabil- 
ity growth rate is infinite. We believe these instabilities 
occur because of coupling in the evolution equations be- 
tween the pure propagating and non-propagating modes 
used in our simple estimate. 
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B. Schwarzschild Spacetime 

We have also studied instabilities in the Einstein 
evolution equations for solutions that are close to 
the Schwarzschild geometry. We use the Painleve- 
Gullstrand [B9l fill, Will form of the Schwarzschild metric: 



ds 2 



-dt z 



dr + y/2M/rdt 



r 2 dn 2 , 



(4.10) 



where dfl 2 represents the standard metric on the unit 
sphere. The dynamical fields that represent this geome- 
try are also quite simple (although not quite as simple as 
Rindlcr). In Cartesian coordinates we have 



9ij 

Ki4 



= 5.. 



D 



kij 

N 



V2Mp% 

0, 

1, 

k 



3y/M/2r 



N k = y/2M/rf 



(4.11) 
(4.12) 
(4.13) 
(4.14) 
(4.15) 



where 8%j is the Euclidean metric (in Cartesian coordi- 
nates), and f k — (x, y, z) jr is the unit vector in the radial 
direction. In this representation of the Schwarzschild ge- 
ometry the fundamental representation of the dynamical 
fields tig = {6ij,Kij,0} includes a non- vanishing extrin- 
sic curvature. Therefore d^T a pu^ will be non-zero for 
this geometry. Since the D^ij components of Uq are still 
zero, it follows that only the Kij components of T a p will 
contribute to V^ a . One can show that only the Kij com- 
ponents of V^aSuQ are non-zero, and that these are given 
by: 

[(1 + 3z)KSfSj - zKg t3 g ab - 9ij K ab ] 6g ab . 

(4.16) 

Thus V^a depends only on the kinematical parameter z. 
Consequently for evolutions near the representation of 
the Schwarzschild geometry considered here, the growth 
rates of the instabilities will depend on only three of the 
nine KST parameters: {C,7, z}. 

We solve the evolution equations for the perturbed 
Schwarzschild geometry using a pseudospectral colloca- 
tion method (see, e.g. ^?], Q for a general review, 
and |2^, [l5| [l9| for more details of the particular imple- 
mentation used here) on a 3D spherical-shell domain ex- 
tending from r = 1.9M to r = 11.9M. Our code utilizes 
the method of lines; the time integration is performed 
using a fourth order Runge-Kutta algorithm. The inner 
boundary lies inside the event horizon; at this bound- 
ary all characteristic curves are directed out of the do- 
main (into the black hole), so no boundary condition 
is required there and none is imposed ("horizon exci- 



sion |£9 go} |5J |5J, |3] ^ [55J 

This reflects the causality condition that the interior of 
a black hole cannot influence the exterior region. At the 
outer boundary we require that all ingoing characteris- 
tic fields be time-independent, but we allow all outgoing 




FIG. 6: Dashed curve illustrates the simple analytical growth 
rate estimate for the Schwarzschild spacetime with k n — 0, 
and the solid curve shows estimate for k n — —r n /M. Points 
are growth rates determined numerically for long-term (many 
light crossing times) evolutions of the 3D linearized equations 
with £ = — 1 and z = — 1/4. 



characteristic fields to prop agate freely. T he ini tial data 
is the exact solution Eqs. ( 4.11 ) through ( 4.15| ), plus a 
perturbation of the form Ae~^ r ~ 7M ^ / 4M added to each 
of the 30 dynamical variables (the Cartesian components 



of gij, K^, and Mkij). The perturbation amplitude A 
for each variable is a randomly chosen number in the in- 
terval (— 10~ 8 ,10~ 8 ). The gauge fields (3 l and Q are not 
perturbed. Because we perturb the Cartesian compo- 
nents of each field by a spherically symmetric function, 
the initial data are not spherically symmetric. Note that 
we have chosen an initial perturbation that violates the 
constraints. 

At the outer boundary, the modes that appear non- 
propagating to a hypersurface orthogonal observer actu- 
ally are incoming with respect to the boundary because 
of the outward-directed radial shift vector N l . Thus the 
projection operator v needed to construct our growth 
rate estimates from Eq. (2.18) is therefore the one that 
annihilates both the incoming and the "non-propagating" 
modes, but leaves the outgoing modes unchanged. We 
have computed the eigenvalues of this appropriately pro- 
jected C a p, and illustrate the largest eigenvalue in Fig. ||. 
These eigenvalues depend on the radial coordinate r in 
this spacetime, and we plot in Fig. ^ the largest value of 
this eigenvalue, which occurs at r = 2M. The graph rep- 
resents (half) this largest eigenvalue as a function of the 
parameter 7 for £ = — 1 and z = — 1/4. In Fig. || we give 
estimates for two choices of the wavevector k a that char- 
acterizes the spatial variation of Su a : the dashed curve 
represents the choice k n = 0, while the solid curve repre- 
sents the choice k n = —r n /M. The points in this graph 
represent the numerically-determined growth rates of the 
instabilities for the linearized equations. We see that the 
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FIG. 7: Same information as Fig. |6| but plotted here using 
a logarithmic scale. This illustrates that while the simple 
analytical growth rate estimate is qualitatively correct, it fails 
to correctly identify the optimum choice of parameters. 



7 = —16. Again, the actual minimum growth rate is 
smaller than the analytical estimate. However, the qual- 
itative shape of the analytical curve is correct, and can 
be used as a guide for choosing parameters to investigate 
with the numerical code. In the present case this guide 
proves extremely useful, for it has allowed us to find a 
parameter choice (z = —0.42, 7 = —16, ( = —1) that 
significantly extends the amount of time the fully 3D non- 
linear code can evolve a single Schwarzschild black hole. 
The very narrow range of parameters where the evolu- 
tion equations have optimal stability, as shown in Fig. |^, 
illustrate why these optimal values were not discovered 
empirically ^6j. The growth rate estimates derived here 
were needed to focus the search onto the relevant part of 
the parameter space. The improvement in 3D nonlinear 
black hole evolutions resulting from these better param- 
eter choices will be discussed in detail in a forthcoming 
paper. 



V. DISCUSSION 
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FIG. 8: Dashed curve illustrates the simple analytical growth 
rate estimate for the Schwarzschild spacetime with k n = 0, 
and the solid curve shows the estimate for k n — —r n /M. 
Points are growth rates determined numerically for long-term 
(many light crossing times) 3D evolutions with £ = — 1 and 
7 = -16. 



agreement with our very simple analytical estimate is 
again quite good. However, Fig. [7] illustrates that while 
the simple analytical estimate is correct qualitatively, it 
fails to correctly identify the best choice of parameters 
to use for long-term numerical evolutions. The minimum 
growth rate is actually smaller (by about a factor of two) 
than what could be achieved by the optimal range of 
parameters that are identified by the simple analytical 
growth rate estimate. 

Fig. H is the same as Fig. [j] except the results are plot- 
ted as a function of the parameter z for £ — — 1 and 



This paper studies the constraint-violating and gauge 
instabilities of the Einstein evolution equations. We de- 
rive an analytical expression for the growth rate of an 
energy norm in Sec. [h) This energy norm measures 
the deviations of a given solution from a background 
constraint-satisfying solution. We show numerically that 
the growth rate of this energy norm is identical to the 
growth rate of constraint violations in solutions of both 
the linearized equations and of the full nonlinear equa- 
tions. Thus we concentrate here on the analysis of the 
evolution of this energy norm. Section III derives the an- 
alytical expressions needed to evaluate this energy norm 
in the 12-parameter family of Einstein evolution equa- 
tions introduced by Kidder, Scheel, and Teukolsky [26). 
This analysis demonstrates that an open subset of the 
KST equations with all physical characteristic speeds is 
symmetric hyperbolic. And perhaps more surprisingly, 
there is a large open subset of these strongly hyperbolic 
equations that is not symmetric hyperbolic. 

We make a considerable effort in this paper to demon- 
strate that the instabilities that we observe are actual 
solutions to the evolution equations and not numerical ar- 
tifacts of the discrete representation of the solution that 
we use. We run many of the cases that we study here 
using both a finite difference and a pseudo-spectral ver- 
sion of the code, and find equivalent results. We run all 
numerical simulations at multiple spatial and temporal 
resolutions to demonstrate that numerical convergence to 
the desired accuracy is achieved. We also confirm that 
the growth rate determined from our analytical integral 
expression is equal (to the desired numerical accuracy) to 
the growth rate measured numerically from the growth 
of the energy norm. This equality is expected only for 
solutions to the actual evolution equations, so this pro- 
vides a non-trivial check that the instabilities we find are 
real solutions to the equations. 
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The analysis presented here also suggests that the in- 
stabilities that we see are endemic to the Einstein evolu- 
tion equations and are not the result of improper bound- 
ary conditions. We impose conditions that suppress 
the incoming components of the dynamical fields at the 
boundaries of the computational domain. These condi- 
tions (sometimes called maximally dissipative J61[) en- 
sure that the energy norm does not grow due to en- 
ergy being inserted into the computational domain across 
the boundaries. Any other boundary conditions that 
might be imposed (including ones that attempt to con- 
trol the influx of constraint violations [p2[ ) could only 
increase the growth rate of our energy norm. Because 
the constraint violations grow at the same rate as the 
energy norm, it seems very unlikely that the constraint 
violations could be eliminated or even significantly re- 
duced simply by changing the boundary conditions. And 
even if the constraint violations could be eliminated 
with better boundary conditions, our analysis shows that 
something — presumably gauge instabilities — would still 
cause the energy norms to grow on the timescales illus- 
trated here. 

The analysis presented here provides additional sup- 
port for the claim that the growth rate of instabilities is 
strongly affected by the choice of representation of the 
Einstein evolution equations. We find that this growth 
rate varies considerably as the parameters in the KST for- 
mulation of the equations are varied. Further, we demon- 
strate here that the functional dependence of the growth 
rate on the parameters strongly depends on the initial 
data that are being evolved. We show that the func- 
tional dependence on the set of parameters is not the 
same for the Schwarzschild geometry as it is for flat space- 
time in Rindler coordinates. This result strongly suggests 
that analyzing the stability of the evolution equations 
for perturbations of flat spacetime in Minkowski coordi- 
nates |H[ ||(| , although useful for screening out particu- 
larly poorly-behaved formulations, is unlikely to succeed 
in identifying the best form of the equations to use for 
binary black hole spacetime evolutions. Rather these re- 
sults suggest that it may be necessary to choose the opti- 
mal form of the evolution equations individually for each 
problem. Estimates of the instability growth rates such 
as those found here (and hopefully more refined estimates 
yet to be discovered) depend only on the fields at a given 
instant of time. It might be possible (or even necessary) 
then to use such estimates to determine and then adjust 
the form of the evolution equations to be optimal as the 
spacetime evolves. 



to use the Cornell code to test our growth rate estimates. 
Some computations were performed on the IA-32 Linux 
cluster at NCSA. This research was supported in part by 
NSF grant PHY-0099568 and NASA grant NAG5-10707. 



APPENDIX 

In this Appendix we give the details of a number of 
equations that define the KST |?6| formulation of the 
Einstein evolution equations. The principal parts of these 
equations can be written in the form of a first order sys- 
tem: 



8 t u a + A ka p d k u 



0. 



(AT) 



Here the dynamical fields are taken to be the set u a — 
{gij , Pij , Mkij} , where gij is the spatial metric, and Pij 
and Mkij are fields that initially will be interpreted as 
the extrinsic curvature Kij and the spatial derivatives 
of the metric D k %j — ^dk9ij respectively. The matrices 
A ka p can be written quite generally for these systems in 
the form 

dtgij 
dtPij 



N n d n9ij , 



be cn 



N 



hi i 



(A.2) 

fHg nb S c l S d 3 +^g nd 5 b {l S c j} 

be cn cd , cd cn cb 

+M33 (i° j) + M4ff (%0 j) 

•/^//'"V" .</,., + ^ 6 g nb g cd g l3 ] d k M bcd , (A.3) 
N n d n M kij - N^^kS^j + v 2 5 n {i 5 b j) 8 c k 
+^g nb 9k(iS c j) + V4g nb gij5 c k 



+^5fi ' 9k(iS n j) + v 6 g g tl $ n k\ d n P bc . (A.4) 

For the fundamental representations of the theory dis- 
cussed in Sec. [II where the dynamical fields are Uq = 
{gij, Kij, Dkij}, the constants ma and ua are related to 
the 5 dynamical KST parameters (7, C, r\, x, cr) by, 



Ml 


= Vl = l, 


(A.5) 


/'2 


= -i-c, 


(A.6) 


/'3 


= -i + c, 


(A.7) 


t'-l 


= 1 + 2(7, 


(A.8) 


M5 


= -Me = -7, 


(A.9) 


v% 


= 0, 


(A.10) 


vz 




(ATI) 




= -ve = -ix- 


(A.12) 
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KST |2q| generalize this fundamental representation of 
the theory by introducing a 7-parameter family of trans- 
formations on the dynamical fields. These transforma- 
tions define the new fields Pij and Mkij in terms of the 
fundamental fields and Dkij via an equation of the 
form 



pu . 



(A.13) 
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The expli cit for m o f this transformation is given in 
Eqs. ( 3.12 ) and ( 3.13 ). Here we give explicit expressions 
for the inverse transformation: 



K l3 =P tJ +zg tJ g ab P ab , (A. 14) 

D kij = [kStfiSj + eS^St +a 9l] g bc S a k 

+bg ij9 ab Sl + c 9k{i 5 a j)9 bc + d 9k{i 6 c j)9 ab ]M abc , 

(A.15) 

where the constants {z, k, a, b, c, d, e} are functions of the 
kincmatical parameters {z, k, a, b, c, d, e}: 



z — 

5a — 

5b = 

5c = 

5d = 

5 e = 

5 a k = 

So = 

5 = 



-z/(l + 3z), (A.16) 
2(3e + 4k)(bc - ad) - [a-b- c + d)e 2 

-2{2de -be - ce + 2ak)k, 
4(2e + k) (ad - be) + 2 (a - c)e 2 

+2(2ae - be + de - 2bk)k, 
4(2e + k){dd - be) + 2 (a - b)e 2 

+2(2ae - ce + de - 2ck)k, 
4(e + 3k)(bc - ad) - 4ae 2 

+4(be + ce - dk)k, 

2e, 

-e- 2fc, 
e 2 - ek-2k 2 , 

5 [I0(6c - ad) + (36 + 3c - a + d + e)e 

-(6a + 26 + 2c+4d + e + 2fe)jfe]. (A.24) 



(A.17) 

(A.18) 

(A.19) 

(A.20) 
(A.21) 
(A.22) 
(A.23) 



The transformation T"^ is invertible as long as 5 ^ 
and z ^ — | . In this generic case we may write 



(A.25) 



and T a 1 T~ l1 p = 5 a p, where 5 a p is the identity matrix 
on the space of fields. 

We have seen in Sec. HI that the general expression 
for the matrices A p is related to its form in the funda- 
mental representation of the equations, Aq u p, by the sim- 



ple transformation law: A ka p = T a ^A^ v T~ Vv p. This 



transformation may also be expressed somewhat more 
concretely as a transformation on the constants \xa and 
va that define the matrices A ka p. The resulting expres- 
sions for these constants after the kinematical transfor- 
mation are: 

M = fc-i(l + C)e, (A.26) 
= i(l-C)e-(l + C)fc, (A.27) 
Ha = (1 + 6<j)6 - (1 - C)fc - |(l-4cr-3C)d 

+i(l+4a + C)e, (A.28) 
m = (1 + 6<r)a + (1 + 2o-)fc - |(l-4(T-3C)c 

-4(1 -C)e, (A.29) 
/i 5 = (1 + 2 7 + 4z + 6jz)(b- |d) + 2o-z(36 + d + e) 



-(7 + 2z + 37i) (fe - |e) - 4Cd, (A.30) 
fi e = (1 + 27 + 4z + 67z)(a - ic) + 2az(3a + fc + c) 

+(7 + 2z + 3jz) (k - Kc, (A.31) 

i/ x = fc, (A.32) 
U3 = e, (A.33) 

^3 = (1 - »7~ hX)d- i{rj + 3x)c- |(?? + 2x)e- |7?fc, 

(A.34) 

z> 4 = (I-77- |x)6- i(r? + 3x)a- irje - (A.35) 
U5 = 4 (2 + 77 + 3x + 6z + 2r?z + 6x^)c 

+4(277 + x + 2z + 4ryz + 2xz)d + 4(77 + 2?7z)fc 

+ i(77 + 2x + 4z + 2ryz + 4xz)e, (A.36) 
vq = § (2 + 7] + 3x + 6z + 2??z + 6xz)a 

+4(2r/ + x + 2z + 4tjz + 2x^)6 

+ i(77 + 2ryz)e + i(x + 2z + 2xz)fc. (A.37) 

These expressions are identical to those in KST [^6| ex- 
cept for ^5 and /X6, which differ from the the KST ex- 
pressions (due to a typographical error in KST) by the 
substitutions c <-> d. 

Finally we wish to give an explicit expression for the 
transformation that relates a fundamental representation 
of the symmetrizer S^p with the general representation 
S a p- We have seen in Sec. [II B that a symmetrizer 
S a p is determined by a set of constants Ba and Ca, 
and similarly the fundamental representation S a p is de- 
termined by constants B A and C A - The general sym- 
metrizer is related to the fundamental by Eq. ( 3.40 ), 
S a p = T- ltJ - a S° u T- 1,J p. This transformation can be 
represented more concretely as relationships between the 
constants Ba and Ca that define S a p and the constants 
B A and C A that define 5°^. These relations are: 



£1 
B, 

Ci 
C 2 

c 3 



- B 2 , 

= (k + 

= (k-ie 

= A 2 C^ 

= V 2 C^ 

= AVC% 



- B 2 Cl 4 

-£ 2 CU 
VB£C\ 



(A.38) 
(A.39) 
(A.40) 
(A.41) 
(A.42) 
(A.43) 
{A£ + BV)C^, (A.44) 



2ABC%, 
2V£Ct 



where 



A = k + 3a + c, 

B = ie + a + 2c, 

V = e + 36 + d, 

£ = k + he + b + 2d. 



(A.45) 
(A.46) 
(A.47) 
(A.48) 



It is straightforward to verify that these Ba and Ca sat 
isfy the positivity conditions needed to guarantee that 
S a p is positive definite, so long as the B A 
satisfy those conditions. 



and C A also 
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additional constraints as well. We illustrate in Sec 
that all the constraints grow at the same rate. 
This projection operator is not unique, however, in the 
cases that we study here there is a obvious simple choice. 
Our Dkij differs from the quantity dkij used by KST by 
a factor of two: dkij = 2Dkij ■ 

We note that this representation in terms of 7 and £ 
encompasses both cases of the representation in terms of 
T) and C given by KST in their Eqs. (2.38) and (2.39) Q. 
While this assumption is not essential, such a sym- 
metrizer must exist if there is any positive definite sym- 
metrizer that depends only on the dynamical fields u a . 
In a neighborhood of flat space only gtj is non-zero so 
any symmetrizer in this neighborhood must depend on 

ga- 

The common belief that strongly hyperbolic systems are 
also symmetric hyperbolic is true only for equations on 
spacetimes with a single spatial dimension. 
We do not reproduce here the rather complicated expres- 
sions for the eigenvalues of C' a p in the Rindler geometry, 
but we would be happy to supply them to anyone who 
has an interest in them. 



